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We demonstrate a scaling method for non-Markovian Monte Carlo wave-function simulations 
used to study open quantum systems weakly coupled to their environments. We derive a scaling 
equation, from which the result for the expectation values of arbitrary operators of interest can be 
calculated, all the quantities in the equation being easily obtainable from the scaled Monte Carlo 
wave-function simulations. In the optimal case, the scaling method can be used, within the weak 
coupling approximation, to reduce the size of the generated Monte Carlo ensemble by several orders 
of magnitude. Thus, the developed method allows faster simulations and makes it possible to solve 
the dynamics of the certain class of non-Markovian systems whose simulation would be otherwise 
too tedious because of the requirement for large computational resources. 
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I. INTRODUCTION 

The description of the dynamics of open quantum sys- 
tems has attracted increasing attention during the last 
few years 0. The major reason for this is the identifi- 
cation of the phenomena of decoherence and dissipation, 
which characterize the dynamics of open quantum sys- 
tems interacting with their surroundings as a main 
obstacle to the realization of quantum computers and 
other quantum devices Q ■ Secondly, recent experiments 
on engineering of environments Q have paved the way 
to new proposals aimed at creating entanglement and 
superpositions of quantum states exploiting decoherence 
and dissipation 

A common approach to the dynamics of open quan- 
tum systems consists in deriving a master equation for 
the reduced density matrix which describes the temporal 
behavior of the open system. The solution for the mas- 
ter equation can then be searched by using analytical or 
simulation methods, or the combination of both. 

This article concentrates on the developing of new 
Monte Carlo simulation methods for non-Markovian 
open quantum systems. The general feature of the 
Monte Carlo methods is the generation of an ensemble 
of stochastic realizations of the state vector trajectories. 
The density matrix and the properties of the system of 
interest are then consequently calculated as an appropri- 
ate average of the generated ensemble. 

Some common variants of the Monte Carlo methods 
for open systems include the Monte Carlo wave-function 
(MCWF) method 0, @, the quantum state diffusion 
(QSD) [3, fToL ITU ], and the non-Markovian wave func- 
tion (NMWF) formulation unravelling the master equa- 
tion in an extended Hilbert space [lEtim. The MCWF 
method has been very successfully used to model the laser 
cooling of atoms. Actually, 3D laser cooling has so far 
been described only by MCWF simulations Q. QSD 
in turn has been found to have a close connection to 
the decoherent histories approach to quantum mechan- 



ics [T3, and NMWF method has been recently applied 
to study the dynamics of quantum Brownian particles 
[l6L Il7j . The various Monte Carlo methods and related 
topics have been reviewed e.g. in Refs. [sl ITsl IT^. l2fj| 

In general, simulating open quantum systems is a chal- 
lenging task. It has been shown earlier that the methods 
mentioned above can solve a wide variety of problems. 
Nevertheless, sometimes there arise situations in which 
the complexity of the studied system or the parameter 
region under study makes the requirement for the com- 
puter resources so large that the solution may become 
impossible in practice, though not in principle. Thus, it 
is important to assess the already existing methods from 
this point of view, and develop new variants to improve 
their applicability. This is the key point of this article. 

Here, we address the Monte Carlo simulation methods 
for the short time-evolution of non-Markovian systems 
which are weakly coupled to their environments. In this 
case, the dynamics of the system may exhibit rich fea- 
tures, whereas the weak coupling may lead to extremely 
small quantum jump probabilities, the consequence being 
unpractically large requirement for the size of the gener- 
ated Monte Carlo ensemble. To overcome this problem, 
we present below a method which in general allows to 
reduce the ensemble size. 

By studying the Hilbert space path integral for the 
propagator of a piecewise deterministic process (PDP) 
[lj , we show that part of the expectation value of an ar- 
bitrary operator A as a function of time t, (A)(t), has 
scaling properties which can be exploited in Monte Carlo 
simulations to speed up the generation of the ensemble, 
in the optimal case by several orders of magnitude. We 
derive a scaling equation, from which the result for (A) (t) 
can be calculated, all the quantities in the equation be- 
ing easily obtainable from the scaled Monte Carlo simu- 
lations. 

We concentrate first on the Lindblad-type non- 
Markovian case which can be solved by the standard 
MCWF method, and then focus on the non-Lindblad- 
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type case which requires the use of the NMWF simula- 
tions in the doubled Hilbert space. 

The paper is structured as follows. Section [H] intro- 
duces the master equation, the corresponding stochas- 
tic Schrodinger equation, and the appropriate simulation 
schemes for the Lindblad- and non-Lindblad-type sys- 
tems. The Hilbert space path integral method is then 
used to calculate the expectation value of an arbitrary 
operator setting the scene for the scaling method which is 
presented in Sec. IIIII Section lTVI shows explicitly how the 
scaling can be implemented and demonstrates the usabil- 
ity of the method, for the example of quantum Brownian 
motion. Finally Sec. presents discussion and conclu- 
sions. 



II. DYNAMICS OF NON-MARKOVIAN 
SYSTEMS 

We describe first in Sec. Ill Al the master equation for 
the Lindblad-type systems and the corresponding stan- 
dard MCWF method. We then continue in Sec III Bl 
with the description of the non-Lindblad-type case with 
the corresponding stochastic Schrodinger equation and 
NMWF unravelling in the doubled Hilbert space. The 
last subsection lll Cl oresents the calculation of the expec- 
tation value of an arbitrary operator A with the Hilbert 
space path integral method which paves the way for the 
scaling procedure. 

We begin by considering master equations obtained 
from the time-convolutionle ss p rojection operator tech- 
nique (TCL) of the form 

^p(t) = A(t)p(t)+p(t)BHt) 

+ Y J C i {t) P {t)D\{t), (1) 

i 

with time-dependent linear operators A(t), B (t) , Ci (t) , 
and Di (t). 

A. Lindblad-type case: master equation and 
MCWF method 

A specific case of the master equation {Q is the one of 
Lindblad-type [UIHHI 

j t p(t) = -i[H S ,p(t)} + J2^(mLiP(t)Lt- 

i ^ 

\L\L lP {t) - \p{t)L\L-\ , (2) 

where Hs is the system Hamiltonian, 7i(i) the time de- 
pendent decay rate to channel i, and Li is the correspond- 
ing Lindblad operator. 

We define this non-Markovian master equation to be 
of Lindblad-type when the time dependent decay coeffi- 
cients ji(t) > for all times t, and non-Lindblad type 



when 7i(£) acquire temporarily negative values during 
the time-evolution |23|. The Lindblad-type case can be 
treated with the standard MCWF method introduced in 
this subsection and the non-Lindblad case with the 
NMWF method in the doubled Hilbert space presented 
in the following subsection 0, 0] . 

The core idea of the standard MCWF method is to 
generate an ensemble of realizations for the state vector 
ip(t) by solving the time dependent Schrodinger equation 

ttMO = H (t)m, (3) 

with the non-Hermitian Hamiltonian H (t) 

H{t)=H s (t)+H DEC (t), (4) 

where Hs(t) is the reduced system's Hamiltonian and the 
non-Hermitian part HjjEc(t) includes the sum over the 
various allowed decay channels i, 

iW(*) = -y£7iW£l£i, (5) 

i 

where the jump operator Li for channel i coincides with 
the Lindblad operator appearing in the master equation 

During a discrete time evolution step of length St the 
norm of the state vector may shrink due to hi dec- The 
amount of shrinking gives the probability of a quantum 
jump to occur during the short interval St. Based on 
a random number one then decides whether a quantum 
jump occurred or not. Before the next time step is taken, 
the state vector of the system is renormalized. If and 
when a jump occurs, one performs a rearrangement of the 
state vector components according to the jump operator 
Li, before renormalization of ip. 

The jump probability corresponding to the decay chan- 
nel i for each of the time-evolution steps St is 

Pi(t) = SbyiWMLlLiW). (6) 

The expectation value of an arbitrary operator A is then 
the ensemble average over the generated realizations 

1 N 

(A)(t) = -J2(^\ A \^)> ( 7 ) 

i=l 

where N is the number of realizations. 



B. Non-Lindblad-type case: Stochastic Schrodinger 
equation and NMWF method in the doubled 
Hilbert space 

The solution of the general master equation can be 
obtained by using the NMWF unravelling in the doubled 
Hilbert space H = 7is © ~Hs where Hs is the Hilbert 
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space of the system p], • The state of the system is 
described by a pair of stochastic state vectors 



'(*) = 



(8) 



such that 8(t) becomes a stochastic process in the dou- 
bled Hilbert space 7i. Denoting the corresponding prob- 
ability density functional by P[6,t], we can define the 
reduced density matrix as 



P (t)= / D6D0*\<j>){i/>\P[O,t] 



(9) 



The time-evolution of 6 (t) can be described as a piece- 
wise deterministic process (PDP) and the corresponding 
stochastic Schrodinger equation reads Q 



d6{t) = -iG(6,t)dt + 
\\0(t)\ 



E 



{l JMOQ-oit)}™*®, (io) 



where the Poisson increments satisfy the equations 
dNi(t) dNj(t) = SijdNi(t), 

E[dNi(t)] = iigT^p dt ' ( U ) 
and the non-linear operator G(9, t) is defined as 

G(0,t)=^F(t) + lY: ]]J f^ r »(«). (12) 

with the time-dependent operators 

' A(t) 



F(t) = 
Ji{t) = 



B(t) 

Ci (t) 
o A (t) 



(13) 



where A(t), B(t), Ci(t), and Di(t) are the operators 
appearing in Eq. IjTjI. 

The deterministic part of the PDP is obtained by solv- 
ing the following differential equation 



i^6(t) = G(9,t), 
and the jumps of the PDP take the form 



'(*) 



'(*) 



Ci(t)ct>{t) 

\Ji{t)6{t)\\ \Di{t)1>{t) 



(14) 



(15) 



Once the ensemble of stochastic realizations has been 
generated, one can then calculate the density matrix of 
the reduced system from Eq. 10. 



C. The Hilbert space path integral for the 
propagator of the PDP and the expectation value of 
arbitrary operators 

For simplicity, we present here the Hilbert space path 
integral for the Lindblad-type case. The derivation of the 
non-Lindblad-type case follows closely the presentation 
below. 

We assume that the initial state of the system is a pure 
state ipo. In this case the propagator T of the PDP (con- 
ditional transition probability) coincides with the proba- 
bility density functional P of the stochastic process 



(16) 



This quantity describes the probability of the system be- 
ing in the state ip at time t when it was in the state ipo 
at some earlier time to. For short time non-Markovian 
evolutions and weak couplings, we assume that the max- 
imum number of jumps per realization is one. Thus, the 
expansion of the propagator T in terms of number of 
jumps contains two terms: deterministic evolution with- 
out jumps and paths with one jump 

T ty, t\^, 0] = T(°) [V>, f |Vo, 0] + T« [V, #o, 0] . (17) 

With the assumptions above, the expectation value of 
an arbitrary operator A at time t can be calculated as Q 

(A)(t) = [ Di>Dr(i>\A\TP)T[i>,t\^ 0> t ] 



DipDip* (i/)\A\ip) 

{r<°> [V, #o, to] + [V, t\ifo, to] } .(18) 

By calculating T^ ' and TW, see Appendix 1X1 we obtain 
for<A)(t) 

(A)(t) = f Di/>Dil>*{il>\A\il>) 



ds^2-fi(s)\\Lig s (ip ) 

i 

x5(ip-g t (ipo)) 

+ J ds J D^D^ J DifaDift 
x8(i>-g t ^ 2 ))^ 7i (s)||L i Vi|| 2 



xS 



Here, terms of the form 5 (%p — g t (ipo)) are the func- 
tional delta-functions and the deterministic evolution of 
tpo according to the non-Hermitian Hamiltonian H is 
given by 



(-i foH(t')dt') xPo 
9t (fa) = ) : ( • (20) 



exp I 



||exp (-if Q H(t>)dt>) Voll' 
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The physical interpretation of Eq. I|19(l is straightforward. 
The expectation value of A is calculated with respect to 
all possible paths of ip with appropriate weights. The 
first term in the curly brackets is the no-jump evolution 
of ip multiplied with the corresponding probability of no- 
jumps. The second term includes the integration over all 
possible jump times and jump routes with the appropri- 
ate transition rates for the one jump realization. 



III. SCALING 

Denoting the expectation value of A with respect to 
the no-jump evolution as 

(A) (t) = (ip = g t (Vo) \A\1> = g t (Vo)>, (21) 
we obtain from Eq. (|19|l 

(A) (i)-(A) (*) = J DipDrmm 

x \ -6 & f dsY,li(s)\\L igs m\ 2 



+ J dsj D^Dipl J Thkrhfe 6 (i> - g t (tp2)) 
Liipi 



xS (tpi - g s (tp )) 



IMill 



(22) 



This equation leads to the first key observation of the 
paper. We notice that (A) (t) - (A) (t) (but not (A) (t) 
alone) is directly proportional to the transition rates of 
the type 



W[V*hM=I>(*)P^i 



\L t ipi 



■ 4>2 



(23) 



In the corresponding Monte Carlo simulations for the 
case we are considering, the required size of the gener- 
ated ensemble is related to the transition rates W since 
the rate defines the number of jumps. In more detail, 
if the total (cumulative) jump probability for the time 
evolution period of interest is P c , we need on average 
to generate 1/P C realizations to produce one realization 
which has a jump. To achieve good statistical accuracy 
we need obviously a large enough number of jumps and 
the minimum condition for the required ensemble size iV 
becomes N 3> 1/P C - 

This leads us to the following observation which can 
be used to optimize the ensemble size of the Monte Carlo 
simulations (within the approximations we use). We can 
artificially increase the number of jumps by scaling up 
the transition rate W by a factor of j3. At the same 
time we must leave the non-Hermitian Hamiltonian H 
unsealed since the ensemble average contribution given 
by realizations with H only (no jumps) appears on the 



l.h.s. of the equation l|22[). In other words, we are not 
allowed to scale the deterministic evolution of the state 
vector (which includes also the rotation of the state vec- 
tor towards the state with the smallest decay rate 7;) 
but only increase the number of jumps by scaling up the 
transition rates by a factor of (3. In the simulation this 
can be done easily multiplying the jump probabilities for 
various decay channels by a same factor (3. An explicit 
example how to do this for both of the cases we are con- 
sidering, Lindblad- and non-Lindblad-type, is shown in 
the next section. 

The question is now how we can calculate from the 
scaled simulations the result we are looking for, namely 
the expectation value for arbitrary operator A as a func- 
tion of time (A)(t). It can be shown, see Appendix [51 
that the final result for {A) (t) starting from Eq. I|22|l can 
be obtained as 



(A)(t) 



1 ~— 8~-0 N ]{A)Q(t) 



+~(*) ! ,,(n. 



(24) 



This equation is the main result of the paper. It shows 
that the ensemble average of the scaled simulations can 
be used to calculate the result for the original problem 
we are interested in. In this equation, P to t (t) is the total 
transition rate (see Appendix [BJ), N is the size of the en- 
semble, Nj{t) the number of jumps in the simulations as 
a function of time, (3 the scaling factor, (A)o(t) the expec- 
tation value with respect the deterministic time evolution 
(see Eq. H211 ). and (A) tot (t) the ensemble average from 
the modified simulations where the scaling has been used 
(see the discussion above). All of the quantities on the 
r.h.s. can be easily calculated in the simulation. Actu- 
ally, from a technical point of view, the only difference 
between the scaled and unsealed simulations is that in 
the former one we have to keep track of the number of 
jumps as a function of time. A task which can be eas- 
ily done in the simulations. We also note that at time 
* = 0, P tot (0) = 0, N d (Q) = 0, (A) tot {0) = (A)o(0) and 
we obtain correctly for time t = 0: (A)(0) = (A)q(0). 

Thus, we can optimize the ensemble size by using the 
following procedure in the Monte Carlo simulations: i) 
Scale up jump probabilities by suitable factor f3. ii) Leave 
decay rates 7$ (t) untouched in the non-Hermitian Hamil- 
tonian H iii) Calculate the result for (A)(t) from Eq. J23J. 

It is worth to emphasize here a common feature of 
Monte Carlo wave-function simulations. The determinis- 
tic evolution caused by the non-Hermitian Hamiltonian 
H changes the relative weights of the occupied states due 
to the different decay rates of the various states. The 
scaling procedure incorporates this rotation by adding to 
the scaled ensemble average result [the second term on 
the r.h.s. of Eq. I|24|l ] contribution from the determinis- 
tic evolution calculated with the appropriate weight (the 
first term). 

Since the reduction in the required ensemble size is di- 
rectly proportional to the used scaling factor /3, the issue 
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is now how large scaling factor we can use to optimize 
the simulations. The scaling method that we have de- 
veloped is valid when there is maximally one jump per 
realization. This condition has to hold also for the scaled 
simulations as well. As soon as the scaling factor is so 
large that realizations with two or more jumps begin to 
occur, additional error (with respect to the normal statis- 
tical error of Monte Carlo simulations) starts to appear. 
In other words, the probability of having two jumps per 
realization has to be much smaller than the one jump 
probability. If the total probability for one jump is P c 
(see the discussion above), the probability for two jumps 
equals P c 2 and the estimate for additional error is sim- 
ply given by P 2 jP c — Pc- Thus we can use the scaling 
factor which increases the jump probabilities e.g. to the 
order of 0.01 introducing a manageable 1% error in ad- 
dition to the normal statistical error of the Monte Carlo 
simulations. 

For the standard Monte Carlo simulations there ex- 
ists a corresponding measurement scheme interpretation 
based on the continuous monitoring of the environment 
of the system. The scaling technique modifies the Monte 
Carlo simulation method in such a way that the measure- 
ment scheme interpretation is lost. 

The scaled simulations correspond to a stochastic 
Schrodinger equation where the deterministic part gener- 
ated by G, see Eq. I|10fl . remains the same but the jump 
part is scaled with 0, i.e. the expectation value of the 
Poisson increment becomes 



E[dNi(t)] = pji{t)\\Liij\\ 2 dt, 



(25) 



Thus the stochastic Schrodinger equation does not have 
a corresponding master equation, and actually does not 
need to have one for the scaling to work. This is because 
we are not looking for two master equations whose results 
are scalable from each other. Rather the key point is to 
modify in a suitable way the equations for the simulations 
in order to make them faster and more efficient. 

Summarizing, we have demonstrated above how the 
scaling works for the Lindblad-type master equation with 
time-dependent but always positive decay coefficients 
7i(t). For this the standard Monte Carlo wave- function 
method can be used 0, H3 • In a similar way, it can be 
shown that the scaling works also for the non-Lindblad- 
type case where 7i(i) may acquire temporarily negative 
values. In this case one needs to use the doubled Hilbert 
space unravelling p|. We show examples of the scaling 
for both of these cases in the next Section. 



IV. EXAMPLES FOR SCALING 

The discussion above shows how it is possible to reduce 
the size of the generated ensemble in the Monte Carlo 
simulations for non-Markovian systems. It is worth not- 
ing that for the Markovian case the scaling is not needed 
because the jump probabilities can be increased trivially 
by increasing e.g. the time step size St in the simulations. 



For the non-Markovian case this does not work because 
the main features of the open system dynamics may be 
given by the time dependence of the decay rates, and St 
has to be kept small compared to the temporal variations 
of the decay coefficients. 

We show below two examples for the scaling. In these 
examples we use the scaling factors (3 = 10 4 and 10 5 while 
the generated ensembles have the sizes of the order of 10 5 . 
In other words, without the scaling, the solution of the 
presented problems would require at least 10 9 ensemble 
members. 

To demonstrate the scaling, we perform the simula- 
tions for the short time non-Markovian dynamics of a 
quantum Brownian particle (damped harmonic oscilla- 
tor) We demonstrate both the Lindblad-type, 
and non-Lindblad type cases. 

The dynamics of a harmonic oscillator linearly coupled 
to a quantized reservoir, modelled as an infinite chain of 
quantum harmonic oscillators, is described, in the secu- 
lar approximatio n, by m eans of the following generalized 
master equation 



dp(t) _ A{t)+T{t) 
dt ~ 2 
A(t)-T(t) 



\2ap(t)a) — a) ap(t) — p(t)a^a\ 



2a) p(t)a — acr p(t) — p(i)aej 1 



(26) 



In the previous equation, a and are the annihilation 
and creation operators, and p{t) the density matrix of 
the system harmonic oscillator. The time dependent co- 
efficients A(t) and T(t) appearing in the master equation 
are known as diffusion and dissipation coefficients, re- 
spectively USUI- 

For an Ohmic reservoir spectral density with Lorentz- 
Drude cut-off, the expression for A(t) is [l(| |25j 



A(t) = 2a 2 kT- 



-(1/r) sm(w i)]}, 



[cos ( Wot) 



(27) 



where the assumption of the high temperature reservoir 
has been used. The dissipation coefficient T(t) can be 
written 



r(t)= 



? 2 

r 2 + 1 



1- 



cos(ui t)-re^ ct sm(uj t) . (28) 



Here, r = uj c /luq is the ratio between the environment 
cut-off frequency lj c and the oscillator frequency a 
is the dimensionless coupling constant, k the Boltzmann 
constant, and T the temperature. When r > 1, the decay 
coefficients A(t) ± T(t) > for all times, and the mas- 
ter equation is of Lindblad type. When r < 1, the de- 
cay coefficients A(t) ± T(t) acquire temporarily negative 
values and the master equation is of non-Lindblad type 
|l(j|. For Lindblad-type master equation, one can apply 
the standard MCWF method (Sec. ITVAl ) where as the 
non-Lindblad-type case requires the application of the 
NMWF method in the doubled Hilbert space (Sec. lTVB|) 
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To demonstrate that the scaling works (in addition to 
the rigorous proof presented above), we compare below 
the results obtained from the simulations to the exact 
analytical results 0, [2(| . 



Lindblad-type master equation and MCWF 
simulations 



For the Lindblad type case we choose parameters 
2a 2 kT/uj c = 1.2 x 1(T 6 , r = 10, a 2 u> /w c = 0.5 x 10~ 8 , 
and the scaling factor f3 = 10 4 . The initial state of the 
system is chosen to be a coherent state |£ = \/2) such 
that, at t = 0, (n) = |£| 2 = 2. We emphasize that the 
present pap er generalizes the scaling method we have 
used in [l|j for initial Fock states to arbitrary system 
Hamiltonians and arbitrary initial states. 

The non-Hermitian part of the Hamiltonian is now 
given by [see Eqs. 10), 0, and pi ] 



x 10" 



H 



if i 



DEC 



{[A(t) - Y(t)]aa) + [A(t) + T(t)}a^a} 



The jump probabilities for each time step St and decay 
channel i are now modified so that the jump probability 
for channel 1 (jump up, absorption of one quantum of 
energy from the environment) is 



(30) 



and for channel 2 (jump down, emission of one quantum 
of energy into the environment) 



P 2 {t) = f3 6t[A(t)+T(t)}(^a 



(31) 



The ensemble average is then calculated in the usual 
Monte Carlo way, as presented in Section III Al and the 
simulation results plugged into Eq. (|24|l to get the final 
result. 

Figure^Jshows the excellent match between the analyt- 
ical curve and the simulations using the scalin g. F or the 
discussion of the analytical solution, see Ref. 161. The 
results confirm once more the validity of the scaling pro- 
cedure and show the short time quadratic non-Markovian 
behavior of the average quantum number (n) = (a^a) of 
the oscillator. Moreover, for the parameters used here, 
the scaling reduces the required ensemble size by a factor 
of 10 4 . The simulation here contains 6 x 10 5 realizations. 



B. Non-Lindblad-type unravelling in the doubled 
Hilbert space 

For the non-Lindblad-type case we choose the following 
parameters 2a 2 kT/uj c = 2.4 x 10~ 6 , r = 0.1, a 2 uj /uj c = 
0.5 x 10~ 8 , and the scaling factor (3 = 10 5 . As initial 
state we choose a superposition of Fock states ip = (|0) + 
|1»/V2. 




FIG. 1: Comparison between analytical (solid line) and 
scaled simulation results (circles) with the bars of the stan- 
dard error for the Lindblad-type case. The figure shows the 
behavior of the expectation value of the quantum number (n) 
as a function of time. The initial state of the system is a 
coherent state |£ = y2). For the parameters used here, the 
scaling reduces the required ensemble size by a factor on the 
order of 10 4 . The simulation here contains 6 x 10 realizations. 



The doubled Hilbert space state vector for the har- 
monic oscillator reads 



6{t) 



4>{t) 



Er=o <l>n(t)\n) 
E„=o *Pn(t)\n) 



(32) 



where </> ra (i) and ip n (t) are the probability amplitudes in 
the Fock state basis. 

By comparing Eq. I|26|l with the master equation JQl, 
the operators A(t) and B(t) have to be chosen as 

A{t) = B(t) = -iuaala- ^ {[A(t) +r(<)]aW 

[A(t)~T(t)}aa^}. (33) 

Accordingly, the operators Ci and D, are 



d(t) = D l {t) = VlA(t) - r(t)j qt, 
C 2 (t) = D 2 (t) = y/\A(t) + T(t)\a 

and the corresponding operators Jj, become 



(34) 



.(35) 



The statistics of the quantum jumps is described by 
the waiting time distribution function F w (t) which rep- 
resents the probability that the next jump occurs within 
the time interval [t, t + r). F w (r), derived from the prop- 
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FIG. 2: Comparison between analytical (solid line) and 
scaled simulation results (circles) with the bars of the stan- 
dard error for the non-Lindblad-type case. The figure shows 
the behavior of the expectation value of the quantum number 
(n) as a function of time. The initial state of the system is the 
superposition of Fock states (|0) + |l))/v2. For the parame- 
ters used here, the scaling reduces the required ensemble size 
atleast by a factor of 10 4 . The simulation contains 6 x 10 5 
realizations. The inset shows (in the same scale as the main 
plot) the poor match between the analytical result (solid line) 
and the simulation result without the scaling (circles) with 
6 x 10 s realizations which is three orders of magnitude larger 
than used in the scaling (see text). 



The results show the excellent match between the exact 
analytical solution and the simulation results using the 
scaling with 6 x 10 5 realizations. Again, the results con- 
firm the validity of the scaling procedure. Moreover, the 
inset shows a very poor match between the non-scaled 
simulations with 6 x 10 s realizations [2]} and justifies the 
claim that the reduction in the ensemble size is at least 
on the order of 10 when the scaling procedure is used. 

The reduction of the ensemble size can be estimated 
also by calculating the maximum jump probability of a 
single realization. In the example considered here, the 
maximum probability is of the order of 1CP 7 , in other 
words on average an ensemble size of 10 7 produces one 
jump event in the unsealed simulations. We estimate 
that one needs several hundreds jumps in the simulations 
to produce accurately the rich dynamical features of the 
heating function displayed in Fig. [3 and consequently the 
requirement for the ensemble size is at least 10 9 without 
the scaling. Thus, the reduction in the ensemble size by 
the scaling method is again found to be at least on the 
order of 10 4 . 

It is interesting to compare the various terms in the 
scaling equation f^jl in the non-Lindblad case. Figure OH 
shows the four terms of the scaling equation Q24J1. One 
can notice that two of the terms practically cancel each 
other and the final result is mostly given by the two terms 
presented in Figs. 03 (a) and (d). 



V. DISCUSSION AND CONCLUSIONS 



erties of the stochastic process, reads 
F w (r) 



1 — exp 



(36) 



where for channel 1 (jump up, the system absorbs a quan- 
tum of energy from the environment) 

*(*) = P '^p^ 1 f> + 1) [\Ut)\ 2 + \Mt)\ 2 ] , 

(37) 

and for channel 2 (jump down, the system emits a quan- 
tum of energy into the environment) 



(38) 

Here, the probabilities are scaled with a factor of (3 ac- 
cording to the scaling scheme presented above. When 
the jump occurs, the choice of the decay channel is made 
according to the factors P\{t) and f^W- The times at 
which the jumps occur are obtained from Eq. 1361) by 
using the method of inversion . 

Figure [21 displays the short time oscillatory non- 
Markovian behavior of the average quantum number (n) . 
This type of behavior is studied in detail in Ref. [lfj. 



We have demonstrated a scaling method for Monte 
Carlo wave-function simulations which can reduce the 
size of the generated ensemble by several orders of mag- 
nitude especially for weakly coupled non-Markovian sys- 
tems. The scaling is based on the notion that once in the 
simulations the jump probabilities are scaled, and the de- 
terministic evolution given by the non-Hcrmitian Hamil- 
tonian left untouched, one can obtain the time evolution 
of the observables of interest from the scaling equation 

The scaling has been used in a restricted form, for a 
specific physical system, in Ref. ^(| . In that case the ini- 
tial state of the system was a Fock state. Here, we present 
a generalized scaling scheme which is able to treat arbi- 
trary initial states of the system and arbitrary Hamil- 
tonians. We emphasize that the scaling method works 
very well for solving the short time dynamics of non- 
Markovian, systems, which bear importance e.g. for the 
decoherence studies for quantum information processing 

n. 

In general, non-Markovian systems, even when they 
are weakly coupled to their environments, can posses rich 
dynamical features despite of the fact that the quantum 
jump probability per stochastic realization is small dur- 
ing the time evolution period of interest (see the examples 
above). This is the key area where the scaling method 
we have presented is useful. The small jump probabili- 
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<%t cq,t 



FIG. 3: Contribution from the various terms of the scaling 
equation fZQ. (a) T A = (A) (t), (b) T B = -Pt ot (t)(A)o(t)//3, 
(c) T C = -{(N~N 3 (t))/((3N)](A) (t), (d) T D = (A) tot (t)/(3. 
The terms has been shifted to start from the same initial 
value for easier comparison. Here A is the number operator 
A — a* a. The final result presented in Fig. 2. is mostly given 
as a sum of the terms displayed in (a) and (d). 



lating the jump probabilities from Eqs. © and l|36f) for 
the time period of interest or by monitoring the number 
of jumps in the simulations. As soon as more than one 
jump per realization in the scaled simulations begin to 
occur, one can estimate the error by calculating the ratio 
between the two-jump and the one-jump probabilities per 
realization. In the examples we have described, we have 
not used very aggressive optimization of the ensemble 
size (though the ensemble size reduction is on the order 
of 10 4 ), and no error has been introduced. This has been 
confirmed by monitoring the jumps in the simulations: 
no two-jump realizations was generated. Thus, the error 
bars displayed in the Figs. an <i HI correspond to the 
usual statistical error (standard deviation) of the Monte 
Carlo ensemble. 

In conclusion, the scaling method has limitations (one 
jump per realization) but it is interesting to note that in 
the region where the method can not be applied (more 
than one jump per realization), it is not needed. This is 
because in this region there already occurs large enough 
number of jumps enhancing the statistical accuracy of 
the simulations. In other words, the problem which the 
scaling solves appears only within the region of validity 
of the method. 



ties due to the weak coupling can lead to the situations 
where the requirement for the size of the generated en- 
semble in the Monte Carlo wave function simulations is 
unconveniently large. In these cases, the scaling method 
can be used to reduce and optimize the generated en- 
semble size for efficient numerical simulation of weakly 
coupled non-Mar kovian systems. 

The scaling method presented here can be used when 
the master equation of the open quantum system can 
be expressed in the general form of Eq. obtained by 
the time-convolutionless projection operator techniques 
(the one-jump restriction still applies, see below). To 
compare our method to the other simulation methods 
for non-Markovian systems one should actually compare 
the validity of the TCL with respect to the methods pre- 
sented e.g. in Refs. [2j| |3(], . Thus, making a rigor- 
ous comparison is an involved task and is left for future 
studies. We initially note here that our method is not 
restricted with respect to the temperature of the envi- 
ronment (while method presented in Ref. |3l| is valid 
for the zero-temperature bath) and is valid, at least in 
principle, to the order used in the TCL expansion of mas- 
ter equation to be unravelled (while method presented in 
Ref. [3(j is post-Mar kovian, i.e. first order correction to 
Markovian dynamics). However, it is worth mentioning 
that the validity of the TCL expansion is crucially related 
to the existence of the TCL generator (see e.g. page 447 
of Ref. 0). 

The scaling method is limited to the cases where there 
is maximally one jump per realization in the generated 
Monte Carlo ensemble. Moreover, it is also important to 
note that the same restriction applies also for the scaled 
simulations. These limits can be easily checked by calcu- 
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APPENDIX A: HILBERT SPACE PATH 
INTEGRAL 

Expanding the exponential waiting time distribution 
F and taking into account the terms corresponding to 
maximum one jump per realization for short times and 
weak couplings, the contribution to the propagator from 
the path without the jumps is |l( 

T< > [</>, # 0! 0] = (1 - F [fa, t]) 5[i>-g t (Vo)] = 

X ~ J ds 5Z^( s )H Li f s (^o)H 2 ^ x 

S^-gt^o)) (Al) 

where 8 (ip — gt (V'o)) is the functional delta- function 
and the deterministic evolution according to the non- 
Hermitian Hamiltonian H is given by 

exp (-i /* H(t')dt>) i> 

9t (iM = ) - t { • (A2) 

|| exp (-if*H(t>)dt>) M 
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where 



where Ptot{t) is the total transition rate 



ih 



H = H s --Y,li{t)Lt L i- 



(A3) 



By using the recursion relation for the propagator T [| 
and neglecting the terms of the order of 7i(i) 2 or higher, 
one can now calculate the contribution of the one jump 
path to the propagator as 

T«h/>,#o,0]= / ds [D^DiPt [d^D^* 



<H^-M^))I>(s)l|£^i|| 2 <5 

i 

8 C0i - 3s Wo)) ■ 



(A4) 



where the transition rate summed over the decay chan- 
nels is 



w[HA] = J2^\\ L ^f 6 



Mi\\ 



i>2 



, (A5) 



The physical interpretation of Eq. (|A4|) is straightfor- 
ward. The integrations sums over the various one jump 
routes and over all the possible jump times. 



APPENDIX B: EXPECTATION VALUE 

In the simulations we scale up the jump probabilities 
by a factor (3, and leave the non-Hermitian Hamiltonian 
as it is [includes also ji(t)], we get the corresponding- 
equation for Eq. ll'L'l) as 



Ptot(t)= I d S V/3 7i (s)||L i5s (^f. (B3) 



Furthermore, we denote by the second term on 
the r.h.s. of Eq. (|B1|) as 



(A>i(t) = / D^DtP* ty\A\tl>) 



ds J D^D^l] D^ 2 Dr 2 
^/5 7l ( s )||^^i|| 2 <5 (ii^ii - V 2 ) S (Vi - g.tya)) 



N 3 (t 



Mlj2(Mt)\A\Mt))/N 3 (t) 7 



(B4) 



where Nj(t) is number of jumps, and N the total number 
of realizations. Here, the second part is the one jump 
contribution to the expectation value, expressed formally, 
and the summation is carried over those realizations that 
have jumped till time t. The corresponding simulation 
presentation (simulation average) is given in the last part. 

Now the ensemble average of all realizations (A) tot (t), 
the quantity which we can easily calculate in the simu- 
lation, is given by as a sum of and 1 jump realization 
contributions 



(A) tot {t) 



N - Nj{t) 
N 



(A) (t) + (B5) 



(3[(A) (t) - (A) (t)} = / DW (1>\A\il>) 



-8(i/> - g t (i> )) J* dsY,01its)\\Li9* 



+ J ds I Di/nDipt I Dip 2 DiP* 2 S (V> - gt(ih)) 



x5 (ipi - g s (ip )) 



\Litpi 



(Bl) 



For scaling to work, we have to be able to extract from the 
simulations the information on the r.h.s. of this equation. 

This can be done as follows. We note the first term on 
the r.h.s. of Eq. (|B1|I as 



(A) (t) = P tot (t)(A) (t) 



(B2) 



Equation I|B1|) , which includes the quantity we are in- 
terested in, can now be written as 



0[{A) (t)-(A) (t)] = -(A) (t) + (A) 1 (t) = 

N — Nj 
— {A + :,) ,,,, 



-Ptot(A) (t) - 



l(A) (t) + (A) tot (t). (B6) 



From this equation we easily obtain the final result for 
the expectation value of arbitrary operator A in compact 
form as 



l-^-i N -.W\{A) {t) 



/3 /3 N 



(B7) 
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